Equilibration of a dissipative quantum oscillator. 



Vinay Ambegaokar 

Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853 

(Dated: May 25, 2006) 

An explicit demonstration is given of a harmonic oscillator in equilibrium approaching the equi- 
librium of a corresponding interacting system by coupling it to a thermal bath consisting of a 
continuum of harmonic oscillators. 

PACS numbers: 03.75.Ss 
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In the theory of interacting many-particle systems, an 
initial condition of thermal equilibrium for the system 
without interactions is often imposed as a convenient de- 
vice, with the implicit assumption that time-evolution 
via the full Hamiltonian will lead to the correct correla- 
tion functions. Such "switching on" of interactions goes 
back a long time 0, > an d is also a part of more recent 
discussions ^ n one form or another, the assump- 

tion is to be found in most theories of transport under 
non-equilibrium conditions Q. However, and despite its 
ubiquity, it is rarely addressed, and sometimes ques- 
tioned [8( . The straightforward and complete demonstra- 
tion within a simple model offered here may therefore be 
of general interest. 

To set the calculation in its general context, consider 
a system consisting of a single particle "sub-system," a 
many-particle "bath," and interactions between the two. 
Let the Hamiltonian for the interacting system be H, 
and let Hq be the Hamiltonian for the complete sys- 
tem with the particle-bath interactions omitted. Now 
consider two reduced density matrices, one describing 
true (time-independent) equilibrium, and a second corre- 
sponding to the time-evolution of the equilibrium of the 
uncoupled system and bath. These may be written as 



p(Q'o,Qo) = TrM){Q' \p th }, and 



(1) 



Po{Qo,Qo,t)=Tr{\Q'i){Q' \ 
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}■ (2) 



Because the couplings in Eq. (3) are quadratic in the co- 
ordinates Qi, where i is or any one of the fcs, the prob- 
lem is formally solved by an orthogonal normal-mode 
transformation 9] to new canonical variables qv,Pv, 



■q v 



1v — E X iv Qi , 



(4) 



with identical connections between Pi and p v . The 
normal-mode frequencies, call them and the trans- 
formation matrices X are obtainable [& 0, 0] from a 
Green function whose inverse is given by 



k 



1 



k\ z 2 



1 , 



(5) 



By putting Eq.(5) on a common denominator, and not- 
ing that the numerator is then the determinant whose 
zeros are the normal mode frequencies, one sees that, as 
a function of the complex variable z, g has the represen- 
tation [lfj 



n*(* 



il(* 2 



(6) 



In words, it has poles on the real axis at the normal-mode 
frequencies, and (interleaved) zeros also on the real axis 
at the bath frequencies. Furthermore, [il flfl Hll | 

1 



Xqu 



In these equations, the primed quantities refer to 
eigenvalues of the sub-system position operator Qq, 
Tr indicates a trace over all states of H, p t h is 
exp{— f3H} /Tr exp{ — [3H} with /3 the reciprocal temper- 
ature, and p\ h is identically defined except that the H of 
Eq.(l) is replaced by Hq. The question being posed is 
whether the quantity defined in (2) becomes (1) for large 
times. 

The model that will be treated here is widely used to 
describe a dissipative quantum oscillator: 

h = \[Pi + n 2 Q 2 Q ] + 1 ]T{ p fe + + KQo/^fc)] 2 }, 

k 

(3) 

with Hq = H(ak — > 0). Here the oscillator is the sub- 
system, and the others the bath. 



1 dg- 
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and 



Xqv 



(7) 



It is convenient to express the projection operator in 
Eqs.(l) and (2) as 



IQoXQol = / dudvf(u,v)e lPoU e lQoV with 



/(«,«) = ±. e -*K* S (Q' -CK-u), 



(8) 



proved by taking matrix elements of both sides, to see 
that I am here asking whether and, if so, how 



Ao(t) = ( e iP °^ u e iQo ^ v ) 



A=(e iPoU e iQoV ). (9) 
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On the left, the operators are in the Heisenberg picture 
corresponding to H and the average is with respect to 
p® h ; on the right, the time- independent average is in p t ^. 
This latter quantity is easily evaluated using the method 
based on the Debye- Waller identity introduced in ref. 
to obtain 

hxA= -±u 2 (P 2 ) - \v 2 {Ql) + $uv, (10) 

For Ao(t), note that the time development due to the 
full H for Pq, Qq is most directly calculated by using 
Eq.(4) to transform to the normal-mode variables p v , q v , 



evolving these according to the quantum mechanical (and 
classical) equations of motion 

Pu(t) = p u (0) coa(w v t) — g , 1/ (Q)w I /sin(o; 1/ t), 

q v {t) = q v (Q) cos(ujJ) +Pu{0)— sm(u v i), (11) 

and then transforming back, using the second part of 
Eq.(4). The Debye- Waller identity may again be invoked, 
because the Schrodinger, i.e. t = 0, Pts and QiS are in- 
dependent oscillator variables for Hq, to obtain 



lnAo(t)= \u^[Y J ^P^^Q 1 )]\-\v 2 ([Y J ^Q l + X:P l )}\ 

i i 



(12) 



Here the Xs are the time-dependent functions 

X?(t) =Y^X 0v X iv coB{u: v t), X°(t) = ^2x 0v X iv —wa.{u v t), and X°(t) = ^X 0l/ X it/ w v wa.(w v t). (13) 

The double sums in Eq.(12) collapse to single sums, again because the oscillators i are independent within Ho, yielding 

i i 

- uvY, [W<J?>o - X?X?(Cfi) a - \iXff - iXfX?]. (14) 

i 

i 



Thus, Eq.(9) is true if the time-dependent coefficients of 
u 2 , v 2 , and uv in Eq.(14) approach the corresponding 
time-independent ones in Eq.(10). 

The sums over v in Eq.(13) are most easily done using 
the method of refs. pi Hoi fll| . In particular, 

XlF{u y ) = 2j>^-zg{z)F{z) 

= 23 / — ug(w - iO + )F(co), (15) 
Jo K 

The contour surrounds the real axis, where the function 
F has been assumed to be regular and zero for lu < 0. In 
the last form the upper limit lu c is the highest normal- 
mode frequency and the integral has been written as the 
imaginary part of the forward portion of the contour. 
Further, using Eq.(7) 



Note that there is no contribution from the poles at 
z = ±f2fc because of the zeros in g(z) at these points, 
evident from Eq.(6). [For an amusing illustration of the 
importance of these zeros, see the final paragraph.] 

Further progress requires a specification of the spec- 
trum of bath oscillators. On general grounds, one would 
expect that to damp out oscillations a continuous spec- 
trum is needed 0] . When this is the case it is meaningful 
to write 

|$>!/fi fe )«Ku;-Q fc ) = J(a;), (17) 
fc 

where J(lu) is a smooth function of its argument in the 
region < u < lu c (w c is the highest bath frequency,) 
and zero elsewhere. Substituting this form into Eq.(5), 
and converting the k sum into an integral, yields 



/dz 
^—2 



g(z) 



l {uj-i0 + ) =uj 2 -Vll-2 



dtu J(ui) 



(z 2 -n 2 k ) 



F(z) 



■> -i diJ g{u-iQ + ) Jr( . 
2a ^ / T%-z0+) 2 -^H- 



(16) 



7T LU (lU — i0+) 2 — LU 2 ' 

(18) 

implying 3g 1 (lu — i0 + ) = —J(uS). 

The time independent part of Eq.(14) may now be ex- 
tracted. [Here I will not examine the remaining time- 



3 



dependent terms, leaving that to further analysis, below, 
of a more specific model.] From Eqs.(13) and (16), 



x c k {t) 



dw g(u> 
— lu— 
ir [lu 



iQ+) 



0+) 2 



«2 



cos(ujt) 



= a k $tg(Slk-iO + )cos(n k t) 



+a k P I — oj- 

7T 



(19) 



In the last line a factor of 2 has been removed and the 
evenness of the integrand used to change the limits of the 
principal value integral (denoted by P.) There is enough 
convergence to set w c -> oo. Now, perform the principal 
value integral by adding one-half of the integral just be- 
low to one half of the integral just above the real lu axis. 
[Note that the scale of 3g(w — i0 + ) is set by J, which is 
not small.] Close the integral in the upper half plane to 
see that 



X^{t) = a k [Kg{n k - ^0 + ) cos(ft fc t) 



-Sg(Q k -iO" 1 



) sin 



(fifc*)]- 



(20) 



where the . . . indicate the time-dependent contributions 
of the singularities of ^sg(uj — i0 + ). The distance of the 
nearest of these singularities from the real axis sets the 
longest time-scale of decaying terms. In a similar way 
one finds 



X^{t) = a k n k [$lg(n k - i0+) sm(n k t) 
+Sg(n k -iO + )cos(n k t)] + 



(21) 



Consider the coefficient of u 2 in Eq.(14). The term 
with i = is not relevant for the present purposes. (See 
below.) Use the thermal averages 

{P 2) =nt(Ql) = ^-coth^, (22) 

to see that the contribution of Eqs. (20,21) for all k to 
the first term on the right hand side of Eq.(14) is the 
time-independent sum 



k 



coth 



|ff(n fc -*o + )| 



(23) 



As the final step, use Eq.(17) to convert the sum over 
k to an integral and exploit — Izp^z" 1 = Qz for any 
complex number to re-express Eq.(23) as 



— 2 W %(c- i 0+)-coth^. (24) 

7T 2 2 



But this is exactly what Eq.(15) gives for the first term 
on the right hand side of Eq.(10), provided that the band 
width of the normal modes, lu Ci is the same as that of the 
bath, lu c . Of course, if this condition were not satisfied, 
the system would not be embedded in the bath and could 
not be expected to equilibrate with it. 



In a similar way, one can show that the other terms in 
Eq.(14) contain time- independent pieces which reproduce 
those in Eq.(10). 

Although there is no time dependence whatsoever in 
these pieces extracted from Eq.(14), it is of course not the 
case that the remainder does not contribute significantly 
at t — 0. In this limit, because of the sinusoidal functions, 
only X? survives. The completeness condition obeyed by 
the orthogonal transformation then requires 



X?(t = 0)=J2 x o»Xiv 



(25) 



making Aq look identical to A, with the important dif- 
ference that the average is in the uncoupled system. 

As mentioned below Eq.(20), the time-decaying parts 
of Ao(t) come from the singularities of ^g(cu — «0 + ), and 
thus depend on the specifics of the bath. It may therefore 
be most instructive to treat a simple model. Consider 
"ohmic" or linear dissipation |7j, in which J(lu) — t\lu for 
< lu < uj c . In this situation Eqs. (17) and (5) give 



z 2 -^-(^)ln(^) 
z 1 - Vl\ + irjz sgn(Sz), |z| < £u c , (26) 



Z + LO c 



so that 



%(w - i0 + ) 



J 2 -fi£)2 +a j2 v 2- 



(27) 



This function has simple poles at the four points ±fi ± 
irj/2, where il 2 = £1q — rj 2 /A. It is these poles that con- 
tribute to the time dependence of the i = terms in 
Eq.(13) and to the parts indicated by the . . .s in Eqs. (20) 
and (21). For an underdamped oscillator, all of them 
thus decay like exp(— rjt/2), showing in this case that the 
damping constant r/ determines the rate at which equi- 
librium is reached. 

This completes the demonstration promised in the ab- 
stract to this paper. 

As an illustration of the importance of the zeros of g(z) 
in the evaluation of sums over normal mode frequencies, 
consider the following proof of the completeness relation 
J2 U X 2 u = 1. Using the method outlined in Eqs. (15) and 
(16) one obtains 



2ai 



dz 



g{z) 



2m {z 2 -fl 2 ) 2 2n k K dz Jz = n «' 



(28) 

The subtracted term is the residue of the simple pole at 
z = Q k which remains after taking into account the zero 
at the same point in g(z). The contour integral can be 
seen to be zero. The second term can be evaluated by 
using an identity expressing the invariance of the deter- 
minant of the dynamical matrix for the coupled oscillator 
system, referred to below Eq.(5), under the orthogonal 



4 



normal mode transformation, 

(z^cid n(- 2 -^)-E «! n = n(- 2 -^)- 

k k k'=ik *> 

(29) 

[Here Og = + J2k a k/^k ^ s tne sc l uare 01 the "bare" 
system frequency in Eq.(3).] Set z — > fi^ and use the 
product representation Eq.(6) for 17(2) to see that the 
term in question is unity. 
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